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SUMMARY 

This  report  describes  a  two  equation  turbulence  model  and  an  alter¬ 
nating  direction  implicit  method  of  solving  the  mean  flow  conservation 
equations  and  turbulence  model  equations,  with  emphasis  on  application  to 
the  calculation  of  the  supersonic  near  wake  behind  a  slender  body.  The 
turbulence  model  has  been  shown  to  adequately  predict  those  low  and  high 
speed  turbulent  flow  experiments  against  which  it  has  been  tested.  The 
numerical  methods  used  have  demonstrated  a  particular  suitability  for  the 
near  wake  problem  because  the  requirement  for  locally  high  spatial  resol¬ 
ution  does  not  require  a  corresponding  reduction  in  temporal  step  size  as 
in  the  case  of  explicit  difference  methods.  Hypersonic  near  wake  results 
are  presented. 
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1.  INTRODUCTION 

The  principal  objective  of  this  program  Is  to  develop  the  capability 
of  performing  calculations  of  turbulent  near  wakes  of  slender  bodies.  Two 
major  subtasks  are  Involved: 

1.  Develop  turbulence  model  equations  applicable  to  compressible 
as  well  as  incompressible  flows 

2.  Develop  finite  difference  methods  for  the  solution  of  the 
model  equations. 

Work  on  the  first  task  has  been  partially  supported  by  the  ROPE  Project 
funded  by  ABMDA  and  ARPA  under  contract  DAHC-G0-71-C-0049  as  well  as  by 
the  present  contract.  An  integral  part  of  the  development  of  the  turbul¬ 
ence  model  equations  Is  comparison  of  calculations  of  simple  flows  with 
experiment.  Such  calculations  are  necessary  to  test  the  validity  of  the 
proposed  modeling  over  a  sufficiently  broad  range  of  conditions  as  to  pro¬ 
vide  some  confidence  In  their  applicability  In  the  near  wake  of  hypersonic 
slender  bodies.  These  comparisons  with  experiment  can  In  most  cases  be 
made  using  subsets  of  the  complete  time-dependent  conservation  and  turbul¬ 
ence  model  equations  proposed  for  use  In  the  near  wake  problem  The  par¬ 
ticular  subset  which  has  been  used  for  this  purpose  is  one  whlci  has  been 
applied  to  the  laminar  steady  calculation  of  interacting  flows/1 }  and 
contains  both  the  complete  invlscld  terms  (in  the  supersede  flow)  and 
boundary-layer-like  viscous  terms.  Since  the  numerical  methods  required 
for  the  implementation  of  the  model  equations  In  this  form  had  already  been 
developed,  this  task  could  be  persued  Independently  of  and  parallel  with 
that  of  developing  suitable  numerical  methods  for  solution  of  the  conplete 
time-dependent  equations.  The  Semi-Annual  Technical  Report(2)  contains  a 
preliminary  description  of  the  model  equations  and  of  the  tests  of  the 
modeling  which  had  been  completed  to  that  date.  Further  tests  have  now 
been  completed  under  the  support  of  the  ROPE  project,  with  a  complete  des¬ 
cription  presented  in  Reference  3.  Agreement  with  experiment  is  now  con¬ 
siderably  better  than  was  shown  in  Reference  2.  The  current  version  of  the 
model  equations  will  be  summarized  in  the  present  report.  Further  testing 
of  the  modeling  will  be  continuing  under  ABMDA  support. 
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The  near  wake  problem,  because  of  the  recirculation  region  immediately 
behind  the  body,  constitutes  a  boundary  alue  problem.  The  proposed  com¬ 
putational  approach,  here,  Is  to  obtain  the  steady  flowfleld  as  the  asymp¬ 
totic  1 1 m t t  of  the  time-dependent  solution  of  the  Kavier-Stokes  equations 
with  turbulence  modeling.  The  numerical  methods  previously  applied  to 
related  problems  have  major  drawbacks  when  applied  to  the  calculation  of 
the  near  wake  flowfleld  at  large  Reynolds  number,  related  primarily  to  the 
fact  that  the  explicit  finite  difference  approximations  used  required,  for 
stability,  chat  the  Courant  number  «c  =  (u  +  c)  be  small'  than  a  limit¬ 
ing  value  which  Is  of  the  order  unity.  (Here, 

u  =  local  velocity 

c  =  local  sound  speed 

At  =  time  increment 

ax  =  local  spatial  meshsize.) 

This  limitation  Is  not  a  serious  one  for  many  applications,  but  proves  to 
be  quite  wasteful  and  inefficient  when  applied  to  the  particular  problem 
of  Interest,  for  two  related  reasons: 

1.  there  are  relatively  limited  flow  regions  which  require  much 
higher  spatial  resolution  than  the  rest  of  the  domain  In  order 
to  avoid  the  generation  of  unacceptably  large  spatial  trunca¬ 
tion  errors  In  these  regions;  for  example,  the  boundary  layer 
and  expansion  fan  In  the  neighborhood  of  the  body  corner. 

Since  the  temporal  step  size  is  severely  limited,  via  the  Courant 
condition,  by  the  small  spatial  meshsizes  required  In  these 
regions  even  when  the  temporal  truncation  error  is  quite  moderate, 
the  explicit  methods  can  be  extremely  inefficient. 

2.  the  approach  to  a  steady  state  occurs  in  two  stages  Involving 
separate  characteristic  times.  For  arbitrary  ini  tic  1  conditions, 
the  flow  rapidly  develops  the  familiar  near-wake  features,  Includ¬ 
ing  the  lip  and  wake  shocks,  expansion  fan,  and  a  re'.lrculatlon 
region.  After  a  time  which  is  of  the  order  of  the  t  anslt  time 
of  a  parti "le  across  the  domain  in  the  outer  high  speed  regions, 
changes  become  much  slower,  and  are  dictated  by  the  development 
of  the  viscous  and  relatively  slow  recirculation  fe^1on.  During 
this  final  slow  approach  to  a  steady  state,  the  e ^ pi  I c i t  methods 
are  still  restricted  In  time  step  by  the  stability  requirements 
associated  with  the  smallest  spatial  mes.<  in  the  expansion  fan, 
which  has  long  since  ceased  to  change  appreciably  with  time  and 
hence  should  not  limit  the  time  step  for  reasons  of  temporal 
truncation  error. 
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For  these  reasons,  an  implicit  finite  difference  method  has  been  devel¬ 
oped  which  is  free  of  these  undesirable  stability  limitations.  This  alter¬ 
nating-direction  implicit  differencing  is  applied  to  the  conservation  form 
of  the  equations  in  order  to  provide  a  better  description  of  the  flow  in  the 
neighborhood  of  shocks  and  other  regions  of  large  gradient.  The  general 
properties  of  these  difference  approximations  applied  to  the  Navier-Stokes 
equations  were  reported  in  the  semi-annual  technical  report,{2)  and  a  des¬ 
cription  of  the  method  with  a  number  of  example  applications  was  presented 

at  the  AIAA  Computational  Fluid  Dynamics  Conference  in  Palm  Sprinqs  ^ 

July  1973. 


These  methods  have  since  been  applied  to  the  problem  of  interest,  the 
hypersonic  high  Reynolds  number  near  wake  of  a  slender  body.  The  remainder 
of  this  report  is  concerned  with  the  description  of  the  application  of  the 
general  methods  described  previously  to  this  particular  prob^m,  along  with 
the  description  of  special  techniques  which  were  found  to  be  required  under 
these  particularly  difficult  conditions,  and  a  presentation  of  the  results 
of  calculations  of  high  Reynolds  number  hypersonic  near  wakes.  The  calcul¬ 
ations  demonstrate  the  particular  suitability  of  the  method  for  this  type 
of  problem  (the  Courant  number  «c  for  these  calculations  was  as  large  as  20 
and  was  never  reduced  below  5).  The  addition  of  a  fl ux-corrected-transport 
(FrT)  stage  t0  the  calculation  is  shown  to  permit  the  description  of  strong 
unresolved  shocks  without  the  severe  numerical  difficulties  normally  assoc¬ 
iated  with  dispersive  effects  which  typically  appear  under  these  conditions. 
This  Is  accomplished  without  resorting  to  artificial  or  numerical  diffusion 
which  can  mask  genuine  viscous  effects  which  are  of  interest,  and  without 
any  increase  in  the  (small)  truncation  errors  associated  with  the  basic 
difference  approximations  as  previously  reported. 
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2.  MODEL  EQUATIONS  KOR  A  TURBULENT  COMPRESSIBLE  FLOW 


The  model  being  used  to  represent  the  non-homogeneous  turbulence  field 
in  the  near  wake  of  a  hypersonic  body  characterizes  turbulence  properties 
at  each  point  by  two  independent  variables  consisting  of  the  turbulent 
kinetic  energv,  e  =  (u'2  +  v'2  +  w'2)/2  and  a  lateral  scale  length  of  the 
turbulence,  i.  To  solve  for  these  variables,  two  partial  differential  equa¬ 
tions  are  written,  each  which  represents  a  bilance  between  the  production, 
dissipation,  diffusion,  and  convection  of  turbulent  energy  and  its  rate-of- 
dissipation.  Because  of  the  closure  problems  always  encountered  in  the 
Reynolds  description  of  turbulence,  the  terms  appearing  in  the  governing 
equations  are  models  of  the  actual  processes  and  are,  therefore,  approximate. 
However,  reasonable  success  has  been  obtained  by  a  number  of  Investigators 
in  the  calculation  of  boundary  layers  and  mixing  layers  and,  consequently 
optimism  exists  for  the  near  wake  problem.  A  more  rigorous  and  detailed 
description  of  the  turbulence  which  considers  its  true  tensorial  character, 
wherein  correlations  between  all  three  components  of  velocity  fluctuation 
as  well  as  orthogonal  integral  scale  lengths  obtained  from  the  various  spec¬ 
tra  are  included,  could  also  be  applied  in  principle.  However,  this  approach 
also  suffers  from  the  uncertainties  of  the  closure  assumptions,  and  the 
additional  complexity  of  many  more  turbulence  variables,  constants  to  be 
determined,  and  equations  to  solve  makes  it  too  unwieldy  for  the  near  wake 
problem.  The  lumped  or  scalar  representation  of  turbulence  by  two  variables 
is  the  simplest  representation  which  contains  sufficient  generality  to  be 
applicable  to  the  problem  of  interest. 

The  mean  flow  equations  are  altered  to  account  for  turbulent  diffusion 
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All  quantities  are  here  non-dimensional i zed  with  reference  length  L, 
velocity  u^,  density  p^,  and  laminar  viscosity  Pressure  Is  normalized 
with  pru£  and  enthalpy  with  ly  E  is  total  energy  CH  -  P/p)  and  H  is  total 
enthalpy  (h  +  +  j-).  Turbulence  is  included  in  the  stress  and  heat 

transfer  terms  through  the  definitions  of  the  transport  properties: 


_  L  .  T 

u  -  U  +  V 


T  1  /2 

where  y  =  pe  '  z  Re 


A,> 


(w) =  (pf)  +(ff) 

2 

where  a  is  normalized  by  ir  and  £  by  L. 

The  governing  equations  for  e  and  £  are  the  turbulent  kinetic  energy 
equation  and  an  equation  for  the  rate-of-disslpation,  e  .  =  c  .  e3/2/£  where 
Cd  is  a  constant.  The  first  equation  has  been  used  for  incompressible  flows  in 
various  modern  investigations,  beginning  with  Bradshaw^  and  is  reasonably 
well  understood  and  accepted.  Less  agreement  exists  with  respect  to  the 
second  turbulence  modeling  equation,  and  the  dissipation  rate  equation 
based  on  the  work  of  Harlow  and  Nakayama^^  and  Jones  and  Launder^3^  was 
chosen  after  comparative  studies  showed  that  this  formulation  agrees  most 
favorably  with  compressible  and  incompressible  boundary  layer  experiments 
(see  References  2  and  3). 


The  turbulence  model  equations  in  the  appropriate  coordinates  are: 
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where  e  =;  CD  e3/2/«,  so  that  uT  -  p  Cq  e2  Re/e 

S]  =  1.0,  S2  =  0.77,  C1  *  1.  C2  =  2.5, 

CD  =  .09,  C]  =  1.41 ,  C2  =  1.89,  PrT  =  0.9 
a  =  local  sound  speed. 

The  laminar  sublayer  dissipation  terms  in  the  turbulence  equations  (des¬ 
cribed  in  References  2  and  3  where  they  are  applied  In  the  compressible 
boundary  layer  calculations  used  to  verify  the  compressibility  modeling 
of  the  diffusion  terms  In  these  equations)  have  not  been  included  here 
since  resolution  of  the  sublayer  on  the  body  base-wall  will  not  be  attempted. 


3.  METHOD  OF  SOLUTION 


The  TACIT  al ternatlng-dl rectlon- Impl left  differencing  method  described 
in  detail  In  References  2  and  4  has  been  modified  slightly  to  Improve  Its 
performance  at  high  Reynolds  number,  when  poor  resolution  In  regions  of 
large  gradient  Introduces  an  associated  large  spatial  truncation  error 
Since  the  method  contains  no  direct  numerical  or  artificial  diffusion  to 
damp  out  such  truncation-error- Induced  oscillatory  disturbances,  they  may 
ecome  large  enough  to  violate  the  requirement  of  positive  static  tempera¬ 
ture  and  terminate  the  calculation.  The  most  direct  way  of  avoiding  these 
problems  Is  to  provide  adequate  spatial  resolution  everywhere.  When  mesh 
spacing  is  uniform,  this  is  not  practical  at  large  Reynolds  number  since 
the  mesh  spacing  required  to  resolve  the  boundary  layers  and  the  expansion 
region  near  the  body  comer  is  so  small  that  the  calculation  becomes  unreason 
ably  costly  In  computer  storage  and  time.  A  second  approach  Is  to  use  a  vari 
ab  e  mesh  spacing,  with  fine  spatial  resolution  used  only  where  required 

This  requires  some  modification  of  the  differencing  presented  In  References 
2  and  4. 

The  variable  mesh  capability  alone  does  not  completely  eliminate  the 
problems  associated  with  local  inadequate  spatial  resolution  at  high  Reynolds 
number,  since  one  does  not  always  know  in  advance  where  all  of  the  problem 
areas  will  occur.  For  instance,  internal  shocks  such  as  the  lip  and  wake 
shocks  in  the  near  wake  cannot  be  precisely  located  a  priori.  For  this 
reason,  a  differencing  technique  introduced  by  Boris  and  Book(9>  which  they 
term  Flux-corrected  Transport  (FCT)  has  been  incorporated  into  TACIT.  This 
technique  consists  basically  of  introducing  an  artificial  diffusion  into 
the  conservation  equations  which  prevents  overshoots  which  would  otherwise 
be  generated  by  large  local  truncation  errors.  Then,  in  a  subsequent  step, 
t  e  effect  of  this  artificial  diffusion  is  removed  identically  everywhere 
except  at  those  few  points  where  this  would  create  a  new  extremum.  At  these 
points,  the  diffusive  flux  used  in  removing  the  effect  of  the  artificial 
diffusion  Is  limited  In  such  a  way  that  the  new  extremum  Is  not  formed, 
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while  still  satisfying  overall  conservation.  The  details  of  the  FCT  algor 
Ithm  Incorporated  Into  a  variable  mesh  version  of  TACIT  will  be  described 
in  Section  3.2. 

3.1  Variable  Spatial  Mesh 


The  terms  whose  difference  approximations  are  affected  by  the  use  of 

g  .  When  these 


a  variable  mesh  are  of  the  general  form  f  and  — 
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terms  are  evaluated  at  meshpoint  1  (corresponding  to  the  spatial  location 
x.|)  the  difference  approximations  used  are: 
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These  general  expressions  replace  those  which  are  presented  in  References 
2  and  4  for  the  special  case  of  uniform  mesh  spacing.  With  these  approxima¬ 
tions,  Inviscid  terms  retain  their  second  order  accuracy  (truncation  error 
of  order  (ax)2)  while  diffusion  terms  are  functionally  second  order  accurate 
provided  Ax  changes  sufficiently  slowly  from  point  to  point. 

3.2  Flux-corrected  Transport  (FCT) 

The  FCT  algorithm  improves  the  behavior  of  the  finite  difference  solu¬ 
tion  in  regions  where  gradients  are  large  and  changing  rapidly.  In  these 
regions,  the  truncation  errors  of  the  difference  approximations  become  large. 
If  the  regions  are  thin  (such  as  in  a  shock  or  thin  shear  layer)  and  the 
difference  approximations  are  conservative,  so  that  overall  conservation 
across  the  thin  region  is  maintained,  these  local  truncation  errors  may  not  be 
of  great  concern.  They  become  important,  however,  if  the  errors  propagate 
away  from  their  region  of  origin  in  the  form  of  wave-like  disturbances  and, 
from  a  practical  standpoint,  these  errors  may  cause  the  calculation  to  fail 
if  they  produce  oscillations  of  amplitude  large  enough  to  violate  a  physical 
limit  to  the  range  of  a  variable  (i.e.,  T  <  0  or  P  <  0).  The  FCT  algorithm 
adds  a  large  artificial  diffusive  terms  to  the  conservation  equations  to 
ensure  that  oscillatory  overshoots  of  this  type  do  not  occur.  The  effect 
of  this  diffusive  term  is  then  removed  in  a  second  calculation  step.  In  the 
absence  of  any  "corrections",  the  net  effect  of  these  two  steps  is  identical 
to  that  of  a  one  step  calculation  in  which  the  artificial  diffusion  term 
is  omitted.  A  "correction"  is  Introduced  during  the  step  which  removes  the 
diffusive  term,  but  only  if  the  effect  of  this  operation  is  to  introduce  a 
new  extremum  In  the  diffused  quantity.  If  a  correction  is  required,  it  is 
imposed  in  conservation  form  (that  is.  If  the  diffusive  flux  leaving  a  volume 
element  is  reduced  as  a  result  of  the  correction,  the  corresponding  flux 
entering  the  adjacent  element  is  correspondingly  reduced  to  maintain  overall 
conservation).  "Corrections"  are  required  only  in  those  local  regions  of 
very  large  spatial  truncation  error  where  the  magnitude  of  the  truncation 
error  Is  not  of  great  concern  since  one  only  expects  correct  overal  1  conser¬ 
vation  across  such  spatially  unresolved  regions.  The  truncation  errors 
associated  with  the  corrections  are  therefore  irrelevant  as  long  as  overall 
conservation  is  maintained.  Within  this  general  framework  describing  the 
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FCT  algorithm,  there  are  many  alternative  detailed  formulations,  of  which 
the  following  has  been  Implemented  in  the  TACIT  code. 

!Hr+-  ;*1/2 *1-1/2) 

(2w2i-i)  r(rJ*''rJ-0 


where 

a 

(7 

■  2i-i)(\+i 

*1+1/2 

~  Sit 

(zi+l 

•Ai) 

*1-1/2  : 

II 

(zi+l  ■ 

•  2i-i)(fli  - 

Ai-l) 

a 

-  ar  / 

{ 

•  rj-i)(Vi 

*j+l/2  ' 

'  2At  ' 

^j+l  ' 

-M 

Vl/2  = 

:  ar  / 

m  i 

VJ+1  “ 

rJ-l)(Aj  • 

AJ-l) 

The  artificial  diffusion  terms  on  the  right  are  differenced  explicitly  or 
implicitly,  consistent  with  the  other  spatial  difference  approximations 

being  used.  The  artificial  diffusion  terms  can  be  considered  as  approxima- 

tion  of  terms  of  the  form  ~  /o,  — \+— —  /Vn  8A  \  ru  .  . 

.  n  at  \U1  9z  /  7  Jr  (rDZ  JF  ■  The  diffusivities  D, 

and  D2  are  functions  of  position  and  are  defined  in  such  a  way  as  to  make 
the  final  difference  approximation  as  simple  as  possible.  The  diffusion 
terms  are  removed  by  means  of  the  following  step: 


a.?.  liVl/2  '  -1/2  )  2(*.M/2  -  *1.1/?)** 

(2i*l  '  Z<-')  rCi*l  -  rJ-l) 

where  I  is  the  value  obtained  when  the  diffusion  term  was  included  in  the 
conservation  equation.  Flux  terms  <p  are  evaluated  at  time  levels  such  that 
the  terms  are  identically  removed  in  the  absence  of  flux  corrections.  The 
corrected  flux  terms  <p  are  evaluated  as  follows: 
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**1+1/2  “  sgn  *1+1/2  max|°* 

(z1+l  "  zi-l) 


Ij.  i  'J1+l 
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(Zi+1  "  z1-l)  /.  \ 

2At 

lM1-l  ■  A1-2  /  sgn  •’i-1/2* 
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W  =  s9"  Vl/2  max  j°- mi"  m  (rj,1  ■  rj-l)(Aj  -  Aj-l)  ^gn 
1 W  -  m  (vi  -  rj.i)(V2  -  Vi)  s9"  tj+i/zj  ] 


*j_l/2  =  s9n  max]°.  min 


24t  (rj+l  ■  rj-l)(Aj-l  •  Aj-2 )  S91 


2ST  (rj+l  -  rj-l)(Vt  '  Aj)sSn  *j.l/2j| 


where  sgn  *  denotes  “the  sign  of  These  tests  constrain  the  change  In  A 
due  to  the  removal  of  the  diffusion  term  in  such  a  way  that  A  can  at  most 
be  made  to  equal  the  values  at  adjacent  meshpoints,  but  cannot  overshoot 
these  values.  Except  In  very  local  regions,  these  tests  do  not  result  In 
corrections,  so  that  almost  everywhere,  J  ■  <p. 
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3.3  Calculation  Sequence 

For  purposes  of  illustration,  we  define  the  following  sequence  of  time 
steps : 


•  »  t 


t  tK+l/2  . 

\  lK+1 


At  the  start  of  the  calculation,  we  assume  thdt  initial  conditions  in  P,  U, 

V,  and  H  are  soeclfied  at  t  =  t, .  Initial  values  of  the  turbulence  properties 
e  and  e  are  assumed  known  at  t  =  t^2.  The  calculation  procedes  as  follows: 


1.  The  mean  flow  equations  are  solved  over  the  interval  t-j  -*•  t2 
using  "odd-step"  differencing  (z  derivatives  differenced  implic¬ 
itly  at  t2»  r  derivatives  differenced  explicitly  at  t1 ,  except 
for  e  and  e  which  are  evaluated  at  t^2) 


2.  The  turbulence  equations  are  solved  over  the  interval  t^^  +  t^2 
using  "odd-step"  differencing  (z  derivatives  differenced  implicitly 
at  t^2,  r  derivatives  differenced  explicitly  at  t^,  except  for 
mean  flow  properties  P,  U,  V,  H  which  are  evaluated  at  t2) 


3.  The  mean  flow  equations  are  solved  over  the  interval  t2  -*•  t3 

using  "even-step"  differencing  (z  derivatives  differenced  explicitly 
at  t2,  r  derivatives  differenced  implicitly  at  t3>  except  for  e 
and  t  which  are  evaluated  at  t^2 


4.  The  turbulence  equations  are  solved  over  the  interval  t 


5/2  7/2 


using  "even-step"  differencing  (z  derivatives  differenced  explicitly 
at  t^,  r  derivatives  differenced  implicitly  at  t^2,  except  for 
mean  flow  properties  which  ore  evaluated  at  t3> 


The  calculation  continues,  using  steps  1-2  if  K  is  odd,  3-4  if  K  is  even. 
The  generation,  dissipation  and  pressure  gradient  terms  in  the  turbulence 
equations,  which  do  not  involve  gradients  of  the  turbulence  variables,  are 
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evaluated  as  average  values  over  the  time  Interval  (half  Implicit,  half 
explicit).  For  it  constant,  the  entire  calculation  sequence  maintains 
second  order  accuracy  In  both  spatial  and  temporal  Increments.  The  split 
t  me  interval  for  mean  flow  and  turbulence  variables  Is  used  for  comput¬ 
ational  efficiency,  since  It  permits  the  sequential  solution  of  (1)  the 
coupled  mean  flow  equations  and  then  (2)  the  coupled  turbulence  equations 
rather  than  requiring  the  simultaneous  solution  of  all  six  coupled  equations. 
The  sequential  solution  Involves  the  solutl.i  of  large  nunfcers  of  coupled 
algebraic  equations  of  diagonal  band  form  with  a  band  width  for  step  (1)  of 
15  elements  and  a  width  for  step  (2)  of  7  elements.  The  simultaneous  solu¬ 
tion  of  al ,  equations  requires  a  band  width  of  22  elements.  Since  the 
Gaussian  elimination  method  used  to  solve  the  equations  requires  computing 
time  proportional  to  m  +  3m  (where  m  is  the  number  of  elements  to  the  left 
or  right  of  the  diagonal,  whichever  is  smaller),  the  sequential  solution 

requires  only  25*  more  computing  time  than  a  laminar  calculation,  while  the 
simultaneous  solution  would  require  85%  more. 

An  additional  advantagt  of  the  sequential  method  of  solution  arises  in 
flow  situations  where  the  turbulence  equations  are  dominated  by  generation 
dissipation  and/or  pressure  gradient  terms.  Under  those  conditions,  the 
mean  flow  equations  require  very  few  iterations  before  the  solution  of  the 
linearized  equations  converges  to  that  of  the  non-linear  finite  difference 
equations.  The  turbulence  equations  may  require  many  more  iterations,  but 
In  the  sequential  method,  each  iteration  of  the  turbulence  equations  takes 

only  25*  of  the  time  of  an  Iteration  of  the  mean  flow  equations.  The  net 
savings  are  very  substantial. 

It  should  be  noted,  here,  that  the  above  situation  (in  which  generation, 
dissipation  and/or  pressure  gradient  terms  dominate)  is  analogous  to  chemi¬ 
cally  reacting  flow  calculations  in  which  a  local  equilibrium  Is  being 
approached,  l.e.,  the  equations  become  stiff  and  hence  unstable  when  solved 
using  explicit  difference  methods  unless  extremely  small  time  Increments  are 
used.  The  Implicit  method  used  here  has  no  difficulties  under  those  condi¬ 
tions,  and,  In  fact,  can  obtain  an  exact  "equilibrium"  solution  (in  which 
convection  and  diffusion  terms,  including  the  time  derivative  term,  are 
identically  zero)  If  reasonably  accurate  Initial  estimator  are  provided. 
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4  HYPERSONIC  NEAR  WAKE  CALCULATIONS 


The  turbulence  modeling  and  structure  of  the  TACIT  code,  as  described 
in  Sections  2  and  3.3,  combine  to  uncouple  the  caVulation  of  the  turbulence 
properties  from  the  calculation  of  the  mean  flow  variables.  That  is,  the 
mean  flow  calculation  for  a  turbulent  flow  differs  from  that  for  a  laminar 
flow  only  in  the  viscosity  law  and  Prandtl  number,  and  the  turbulence  con¬ 
tribution  to  total  viscosity  is  constant  over  any  given  mean  flow  calcula¬ 
tion  time  interval.  This  permits  the  turbulent  mean  flow  calculation  pro¬ 
cedure  to  be  completely  checked  out  independently  of  the  turbulence  model 
equations  by  simply  specifying  a  laminar  viscosity  law. 

4.1  Laminar  nypersonic  Near  Wake 

The  test  case  used  for  preliminary  calculations  (using  the  uniform 
spatial  mesh  version  of  the  code)  was  an  8  degree  sharp  cone  with  a  shoulder 
rounded  to  a  radius  r/Dbgse  =  1/12,  =  21  ,  Re^  =  2.5  x  105,  with  a  cold 

wall  =  .02).  The  viscosity  law  used  was  u/u^  =  T/T^,,  with  Pr  =  1. 

The  flow  around  the  corner,  from  an  Initial  8  degree  wall  inclination  to 
the  point  where  the  wall  inclination  is  -29  degrees,  was  calculated  using 
the  steady  state  code  of  Reference  1. 

This  calculation  could  not  be  made  to  reach  a  steady  state  without 
encountering  a  state  in  which  T  <  0  In  a  local  region  associated  with  the 
wake  of  the  outer  edge  of  the  boundary  layer.  The  problem  is  associated 
with  the  large  spatial  truncation  errors  generated  in  this  region  of  large 
and  poorly  resolved  gradients.  The  initial  boundary  layer  profiles  at  the 
body  corner  were  then  rescaled  to  make  the  boundary  layer  thicker  by  a 
factor  of  five  to  provide  a  less  severe  test  case.  When  problems  still 
persisted,  the  program  was  rewritten  to  provide  a  variable  spatial  mesh. 
Using  the  mesh  shown  in  Figure  1,  behavior  was  considerably  improved,  but 
again  T  <  0  was  eventually  encountered  in  the  wake  of  the  outer  edge  of 
the  boundary  layer  as  it  passes  through  the  corner  expansion  fan.  The 
Incorporation  of  the  FCT  algorithm  finally  eliminated  this  problem  and 
a  smooth  approach  to  steady  state  was  obtained. 


The  boundary  conditions  for  this  calculation  were: 


a)  left-hand  boundary 

1 )  on  base  wall : 

2)  above  corner: 

b)  upper  boundary: 

c)  right-hand  boundary: 

d)  lower  (axis)  boundary: 


u=v  =  C,  H=H  ,  P  from  continuity 
w 

u ,  v,  H,  r  specified 

outflow  (u,  v,  H,  P  extrapolated) 

outflow 


in  =  if.  =  M  =  o  v  =  o 
3r  3r  3r  ’  v  u< 


A  map  of  the  computed  static  pressure  is  presented  in  Figure  2.  At 
the  outer  (large  r/D)  edge  of  the  domain,  the  expansion  fan  can  be  seen  to 
have  an  anomalous  structure.  This  is  due  to  an  approximation  used  in  spec¬ 
ifying  the  inflow  boundary  conditions  above  the  body  shoulder;  i.e.,  the 
inviscid  flow  outside  the  boundary  layer  was  taken  to  be  a  uniform  flow 
corresponding  to  wall  conditions  for  a  Taylor-Maccoll  inviscid  cone  flow- 
field.  Since  the  Mach  angle  is  smaller  than  the  flow  inclination  angle 
in  this  high  Mach  number  flow  region,  the  disturbance  is  carried  upward 
and  actually  leaves  the  domain  through  the  upper  boundary.  Unless  the 
calculation  is  carried  much  further  downstream,  where  the  outer  streamlines 
are  turned  Inwards  by  the  expansion  fan  and  the  influence  of  this  outer 
region  finally  becomes  felt  near  the  axis,  the  discrepancies  in  the  outer 
region  are  of  no  consequence. 

With  the  thick  initial  boundary  layer  used  in  this  test  case,  the  lip 
shock  does  not  appear  as  a  separately  identifiable  structure.  The  lip-wake 
shock  has  its  origin  in  the  pressure  gradients  generated  in  the  recircula¬ 
tion  region,  and  is  further  strengthened  by  the  increasing  pressure  in  the 
near  wake,  to  a  point  about  0.5  X/D  downstream  of  the  base  wall.  At  this 
point  of  maximum  pressure,  P/P^  «  3  at  the  centerline.  The  favorable  pres¬ 
sure  gradient  downstream  of  this  point  rapidly  weakens  the  wake  shock. 

The  corresponding  static  temperature  map  is  shown  in  Figure  3.  Since 
the  body  wall  Is  cold,  the  boundary  layer  temperature  has  a  maximum  value 
at  a  point  well  removed  from  the  wall.  The  decrease  in  temperature  of  the 
wake  of  this  part  of  the  boundary  layer  as  it  passes  through  the  corner 


expansion  fan  Is  evident  In  the  figure.  The  region  of  peak  temperature 
reappears  after  crossing  tne  wake  shock,  but  the  peak  is  slowly  eroded  by 
diffusion  to  cooler  surrounding  gas  at  the  same  time  the  flow  nearer  the 
axis  Is  being  heated  by  dissipation.  At  the  downstream  end  of  the  domain, 
the  static  temperature  has  become  almost  uniform  below  the  shock. 

A  corresponding  map  of  flow  speed  q/u^  is  shown  in  Figure  4.  Separation 
occurs  on  the  base  wall  slightly  below  the  "corner".  Remembering  that  a 
steady  code  was  used  to  compute  the  flow  around  the  rounded  body  shoulder 
to  a  point  where  the  wall  inclination  is  -29  degrees  (the  estimated  ser-n- 
tion  point),  this  indicates  that  this  esti  ,ced  inclination  angle  was  Luj 
small  by  several  degrees.  The  velocities  in  the  recirculation  region  are 
quite  small,  due  to  the  combination  of  large  Re  and  thick  boundary  layer. 

The  thin  boundary  layer  on  the  base  wall  is  evident  from  the  .01  velocity 
contour. 

Figure  5  shows  contours  of  total  enthalpy.  Again,  the  total  enthalpy 
in  the  recirculation  re  ion  remains  quite  small  because  of  the  combination 
of  slow  diffusion  (large  Re)  and  small  gradients  (thick  boundary  layer). 

The  convergence  of  total  enthalpy  contours  to  the  axis  downstream  of  the 
wake  stagnation  point  is  due  to  a  combination  of  the  necking  down  or  the 
streamlines  due  to  the  rapidly  increasing  axial  mass  flux  near  the  axis, 
and  the  effects  of  diffusion  and  dissipation. 

4.2  Turbulent  Near  Wake 

A  calculation  has  been  initiated  unde>*  turbulent  flow  conditions,  cor¬ 
responding  to  an  experimental  investigation  by  Gran^^^.  The  body  is  a 
sharp  three  degree  half  angle  adiabatic  cone  with  a  rounded  shoulder 

(r/^cone  -  ^2),  -  7.56,  Re^p  =  1.7  x  10  .  The  inflow  boundary  conditions 

are  obtained,  as  in  the  preceding  laminar  case,  by  performing  a  steady  flow 
calculation  expanding  the  boundary  layer  around  the  comer  to  a  wall  inclin¬ 
ation  of  -28°.  The  initial  conditions  for  this  calculation  are  obtained  by 
starting  with  estimated  profiles  approximately  two  base  diameters  upstream 
of  the  shoulder  and  computing  the  development  of  mean  flow  and  turbulence 
properties  along  the  body.  Comparison  wi**  the  measured  velocity  profile 
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at  the  corner  shows  acceptable  agreement.  The  Initial  values  of  the 
deDendent  variables  are  specified  arbitrarily,  with  the  objective  of  main¬ 
taining  cotfoatlbll Uy  with  the  boundary  conditions  throughout  the  calcul¬ 
ation,  l.e.,  the  flow  must  never  reverse  and  come  In  through  an  outflow 
boundary  and  should  preferably  remain  always  supersonic  with  respect  to 
this  boundary.  Initial  turbulence  properties  are  specified  to  make  the 
turbulence  energy  and  the  turbulent  viscosity  small  (e  ^  10‘6,  ^  0.5). 

The  calculation  lias  not  yet  been  completed,  but  has  proceeded  far  enough 

(tu»/D  *  °’4)  t0  observe  start  of  the  recirculation  region  as  a  separa¬ 
tion  bubble  on  the  base  wall  just  below  the  shoulder.  The  wake  of  the 

boundary  ayer  extends,  at  this  point  in  the  calculation,  about  0.4  diameters 
downstream  of  the  shoulder. 

The  numerical  methods  outlined  in  this  report  perform  well,  with 
three  ite  -ns  taken  to  insure  convergence  and  the  results  demonstrating 
good  convergence  of  the  rapidly  changing  turbulence  variables  after  two 
iterations. 


I 

I 


23 


5.  CONCLUSIONS 


The  objective  of  this  program  has  been  to  develop  the  methodology  for 
performing  calculations  of  the  near  wake  behind  a  slender  body  in  super¬ 
sonic  flow.  The  major  technical  problems  are  the  development  of  turbulence 
model  equations  which  adequately  describe  a  supersonic  turbulent  flow,  and 
the  development  of  numerical  methods  for  solving  these  equations.  A  two 
equation  moael  of  turbulence  has  been  developed  v/hich  ha?  been  demonstrated 
to  adequately  predict  those  low  and  high  speed  turbulent  flow  experiments 
against  which  it  has  been  tested.  Further  testing  of  the  model  Is  continu¬ 
ing  as  part  of  a  separate  program. 

The  calculation  of  the  two-dimensional  flowfleld  described  by  these 
two  turbulence  model  equations  and  the  four  mean  flow  conservation  equa¬ 
tions  has  been  posed,  in  the  present  study,  as  a  time-dependent  problem. 

The  alternating-direction  implicit  methods  developed  have  been  demonstrated 
on  a  variety  of  problems  of  varying  complexity  ranging  from  one-dimensional 
shocks  to  the  final  problem  of  interest.  The  method  has  been  shown  to  be 
especially  suited  to  this  problem  because  the  temporal  step  size  is  not 
limited  by  the  Courant  criterion  of  explicit  methods.  Specifically,  this 
permits  fine  spatial  resolution  locally  without  a  correspondingly  small 
time  step  being  required.  The  Implicit  nature  of  the  calculation  also 
permits  the  description  of  a  local  equilibrium  state  of  the  turbulence 
properties,  which  can  become  difficult  and  time  consuming  when  done  explic¬ 
itly  because  the  equations  can  then  display  an  instability  similar  to  that 
of  the  stiff  condition  encountered  In  chemically  reacting  flows. 

Results  have  been  presented  for  a  hypersonic  laminar  near  wake  and  a 
turbulent  calculation  is  well  under  way.  The  detailed  results  of  this 
final  calculation  will  be  presented  in  an  addendum  to  this  report,  to  be 
distributed  at  a  later  time. 


♦this  work  has  been  partially  supported  by  ABMDA^ 
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